Packages
knitr::opts_chunk$set(
echo = TRUE,
message = FALSE,
warning = FALSE
)
`%nin%` = Negate(`%in%`)
if (!require("tidyverse")) install.packages("tidyverse")
library("tidyverse")# for dplyr workflow of calculations
if (!require("lmerTest")) install.packages("lmerTest")
library("lmerTest") # LMMs
if (!require("ggplot2")) install.packages("ggplot2")
library("ggplot2") # plots
if (!require("RColorBrewer")) install.packages("RColorBrewer")
library("RColorBrewer") # pretty colors
if (!require("ggpubr")) install.packages("ggpubr")
library("ggpubr") # for multi-ggplot figs
if (!require("Hmisc")) install.packages("Hmisc")
library("Hmisc") # correlation matrix
if (!require("PerformanceAnalytics")) install.packages("PerformanceAnalytics")
library("PerformanceAnalytics") # correlation matrix
if (!require("emmeans")) install.packages("emmeans")
library("emmeans") # marginal means & confints
if (!require("car")) install.packages("car")
library("car") # Anova
if (!require("rmdformats")) install.packages("rmdformats")
library("rmdformats") # clean html R markdown formatSet Colors & Themes
savs_ggplot_theme <- theme_classic() +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text.align = 0,
legend.position = "bottom"
)
savs_ggplot_theme2 <- theme_classic() +
theme(text = element_text(color = "black",
family = "sans",
size = 10),
axis.text = element_text(color = "black",
family = "sans",
size = 10),
legend.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text.align = 0,
legend.position = "bottom"
)
my_brew_cols_4 <- brewer.pal(11, "RdYlBu")[c(5,4,10,11)] # this one!Background
This data was collected on Blunt-nosed Leopard Lizards (Gambelia sila) during their active season in 2021. In this document, we will compute metrics of thermoregulation and behavior, then assess whether differences in those variables could be related to differences in hydric physiology.
All of the data wrangling and preparation for the thermal and microhabitat use data is done within this document, but is done in the ‘wrangling_’ documents for the hydric data.
First, we compute maximum temperatures and thermoregulatory accuracy (thermal metrics) and microhabitat use for each lizard at different time intervals, then we pair that with our hydric physiology data, and look for any relationships.
Load & Wrangle Data
Lizard Data
This is the same dataframe used in ‘analysis_1’. For info on data preparation, see the ‘wrangling_’ files, and for an explanation of each of the variable names, see that section in ‘analysis_1’.
# measurement dates
apr_may_dates <- as.Date(c("2021-04-23", "2021-04-24", "2021-04-25"))
may_dates <- as.Date(c("2021-05-07", "2021-05-08", "2021-05-09"))
liz_dat <- read_rds("./data/G_sila_clean_full_dat.RDS") %>%
# we only use data for radio collared lizards in this analysis
dplyr::filter(complete.cases(radio_collar_serial)) %>%
# add date chunks to help look at values over time
mutate(date_chunks = factor(case_when(capture_date %in% apr_may_dates ~ "Apr 26 - May 6",
capture_date %in% may_dates ~ "May 10 - 20")),
radio_collar_serial = factor(radio_collar_serial))
summary(liz_dat)## capture_date_time capture_date radio_collar_serial
## Min. :2021-04-23 10:39:00.00 Min. :2021-04-23 229-069: 3
## 1st Qu.:2021-04-23 14:54:30.00 1st Qu.:2021-04-23 245-742: 3
## Median :2021-04-24 10:50:30.00 Median :2021-04-25 245-744: 3
## Mean :2021-04-28 14:03:08.27 Mean :2021-05-12 245-745: 3
## 3rd Qu.:2021-05-07 11:41:30.00 3rd Qu.:2021-05-08 245-748: 3
## Max. :2021-05-08 14:19:00.00 Max. :2021-07-14 252-883: 3
## NA's :18 (Other):58
## PIT_tag_ID individual_ID sex_M_F mass_g SVL_mm
## Length:76 F-12 : 3 Female:32 Min. :25.00 Min. : 85.0
## Class :character M-09 : 3 Male :44 1st Qu.:33.70 1st Qu.: 96.0
## Mode :character M-10 : 3 Median :38.00 Median :100.0
## M-11 : 3 Mean :39.20 Mean :101.7
## M-19 : 3 3rd Qu.:43.45 3rd Qu.:108.0
## M-20 : 3 Max. :58.80 Max. :122.0
## (Other):58
## tmt tmt_mass_change_g capture_to_msmt hematocrit_percent
## Water Tmt:38 Min. :-3.0000 Min. : 11.0 Min. :17.00
## Sham Tmt :38 1st Qu.:-1.0000 1st Qu.: 62.5 1st Qu.:30.00
## No Tmt : 0 Median : 0.0000 Median : 152.0 Median :34.00
## Mean :-0.2262 Mean : 401.1 Mean :34.35
## 3rd Qu.: 0.0000 3rd Qu.: 265.2 3rd Qu.:38.50
## Max. : 3.0000 Max. :1665.0 Max. :58.00
## NA's :15 NA's :18 NA's :1
## cloacal_temp_C CEWL_g_m2h msmt_temp_C msmt_RH_percent
## Min. :24.00 Min. : 0.650 Min. :19.07 Min. :12.20
## 1st Qu.:30.00 1st Qu.: 7.989 1st Qu.:28.02 1st Qu.:14.79
## Median :32.00 Median : 9.620 Median :29.31 Median :20.38
## Mean :31.92 Mean : 9.854 Mean :29.18 Mean :22.70
## 3rd Qu.:34.00 3rd Qu.:12.324 3rd Qu.:31.74 3rd Qu.:27.49
## Max. :36.00 Max. :21.673 Max. :33.55 Max. :40.22
##
## e_s_kPa msmt_VPD_kPa osmolality_mmol_kg_mean month week
## Min. :2.205 Min. :1.415 Min. :309.3 April:39 1 :39
## 1st Qu.:3.784 1st Qu.:2.487 1st Qu.:345.0 May :23 3 :23
## Median :4.076 Median :3.229 Median :361.0 July :14 12:14
## Mean :4.100 Mean :3.212 Mean :363.6
## 3rd Qu.:4.685 3rd Qu.:3.982 3rd Qu.:377.9
## Max. :5.188 Max. :4.319 Max. :436.5
## NA's :1
## week_num ultrasound_date gravid_Y_N number_eggs
## Min. : 1.000 Min. :2021-04-24 Gravid :16 Min. :3.000
## 1st Qu.: 1.000 1st Qu.:2021-04-24 Not Gravid:12 1st Qu.:3.750
## Median : 1.000 Median :2021-04-24 NA's :48 Median :4.000
## Mean : 3.632 Mean :2021-04-29 Mean :3.812
## 3rd Qu.: 3.000 3rd Qu.:2021-05-08 3rd Qu.:4.000
## Max. :12.000 Max. :2021-05-08 Max. :5.000
## NA's :48 NA's :60
## largest_egg_diameter_cm largest_egg_circumference_cm stage
## Min. :0.5700 Min. :1.710 small round: 1
## 1st Qu.:0.8425 1st Qu.:2.385 large round: 9
## Median :0.9950 Median :2.885 soft : 2
## Mean :1.1131 Mean :3.148 soft oblong: 1
## 3rd Qu.:1.4200 3rd Qu.:3.985 firm oblong: 2
## Max. :1.8000 Max. :4.740 hard oblong: 1
## NA's :60 NA's :60 NA's :60
## dev_point SMI date_chunks
## Min. :1.00 Min. :23.67 Apr 26 - May 6:39
## 1st Qu.:2.00 1st Qu.:32.29 May 10 - 20 :23
## Median :2.00 Median :34.40 NA's :14
## Mean :2.75 Mean :34.77
## 3rd Qu.:3.25 3rd Qu.:37.43
## Max. :5.00 Max. :45.63
## NA's :60# subsets by date
liz_dat_apr_may_only <- liz_dat %>%
dplyr::filter(complete.cases(date_chunks))
liz_dat_july_only <- liz_dat %>%
dplyr::filter(month == "July")
# ID info only
liz_IDs_only <- liz_dat %>%
group_by(radio_collar_serial, individual_ID, sex_M_F) %>%
summarise(n = n()) %>%
dplyr::select(-n)Validation
This is a dataframe I can use to make sure we do not have data for these lizards after these dates, because they were either found dead or their dropped collars were found.
check_gone <- data.frame(Serial = factor(c("229-072", # found dead June 7
"229-060", # collar found May 18
"229-074", # collar found April 28
"229-070", # became -a on May 8
"245-762", # became -a on May 8
"229-073", # dropped sometime before May msmts
"245-751" # died May 7/8
)),
date_def_gone = as.Date(c("2021-06-07", "2021-05-18",
"2021-04-28", "2021-05-07",
"2021-05-07", "2021-05-07",
"2021-05-07"),
format = "%Y-%m-%d"))
check_gone## Serial date_def_gone
## 1 229-072 2021-06-07
## 2 229-060 2021-05-18
## 3 229-074 2021-04-28
## 4 229-070 2021-05-07
## 5 245-762 2021-05-07
## 6 229-073 2021-05-07
## 7 245-751 2021-05-07Tsurf Wrangling
We measured the surface body temperature of lizards throughout the study using Holohil temperature-sensitive transmitters, which is the data we load in this section.
variables in this dataframe:
- DateTime = date and time the temperature was logged
- Serial = the radio-collar serial number for that observation (to be used to link to specific individuals)
- Freq = what radio frequency the collar was on
- pp = the pulse period recorded by the radio-collar receiver (pulse period is a proxy for temperature)
- ppTempC = the pp value converted into temperature in degrees Celsius using third order polynomial functions estimated from the reference transformations provided by Holohil
temps <- read.csv("./data/telemetry/tower_temps_all.csv") %>%
# properly format data classes
mutate(date_only = as.Date(substr(DateTime, 1, 10), format = "%Y-%m-%d"),
hour_of_day = as.numeric(substr(DateTime, 12, 13)),
DateTime = as.POSIXct(DateTime, format = "%Y-%m-%d %H:%M:%S"),
Serial = (paste(substr(Serial, 1, 3),
substr(Serial, 4, 7), sep = "-")),
Freq = round(Freq, 2))
# look at temps dataframe
summary(temps)## DateTime Serial Freq
## Min. :2021-04-25 13:11:35.00 Length:107452 Min. :150.0
## 1st Qu.:2021-05-02 18:18:19.75 Class :character 1st Qu.:150.4
## Median :2021-05-11 07:58:47.00 Mode :character Median :151.1
## Mean :2021-05-15 02:34:25.98 Mean :151.0
## 3rd Qu.:2021-05-17 22:48:35.50 3rd Qu.:151.5
## Max. :2021-07-13 13:11:42.00 Max. :152.0
## ppTempC pp date_only hour_of_day
## Min. :-23.99 Min. :1001 Min. :2021-04-25 Min. : 0.00
## 1st Qu.: 31.80 1st Qu.:1565 1st Qu.:2021-05-02 1st Qu.:10.00
## Median : 36.26 Median :1703 Median :2021-05-11 Median :13.00
## Mean : 35.20 Mean :1734 Mean :2021-05-14 Mean :12.69
## 3rd Qu.: 39.06 3rd Qu.:1885 3rd Qu.:2021-05-17 3rd Qu.:16.00
## Max. : 76.93 Max. :2799 Max. :2021-07-13 Max. :23.00## [1] "229-059" "229-060" "229-063" "229-065" "229-066" "229-069" "229-070"
## [8] "229-072" "229-073" "229-074" "229-077" "229-078" "229-079" "229-080"
## [15] "245-741" "245-742" "245-745" "245-748" "245-751" "245-752" "245-753"
## [22] "245-754" "245-756" "245-757" "245-759" "245-760" "245-761" "245-762"
## [29] "252-881" "252-882" "252-883" "252-884" "252-885" "252-886" "252-887"Check whether we have observations on lizards that were not alive or were not actually wearing their radio-collars because they dropped them.
temp_check <- temps %>%
left_join(check_gone, by = 'Serial') %>%
dplyr::filter(date_only > date_def_gone) %>%
arrange(Serial) %>%
dplyr::select(Serial, date_only, date_def_gone)
unique(temp_check$Serial)## [1] "229-060" "229-070" "229-072" "229-073" "229-074" "245-751" "245-762"Great… we have data for all 7 individuals that we need to clean up for
temps2 <- temps %>%
# delete the obs for all Serials but 229-070 & 245-762
dplyr::filter(!(Serial == "229-072" & date_only >= "2021-06-07")) %>%
dplyr::filter(!(Serial == "229-060" & date_only >= "2021-05-18")) %>%
dplyr::filter(!(Serial == "229-074" & date_only >= "2021-04-28")) %>%
dplyr::filter(!(Serial == "229-073" & date_only >= "2021-05-07")) %>%
dplyr::filter(!(Serial == "245-751" & date_only >= "2021-05-07")) %>%
# change the serial and ID of F-08 and M-03 to be "-a" bc collar was switched
# for F-08a
mutate(Serial = replace(Serial, # what I want to change
Serial == "229-070" & # if this
date_only > "2021-05-07", # and this
"229-070a"), # change to this
# for M-03a
Serial = replace(Serial, # what I want to change
Serial == "245-762" & # if this
date_only > "2021-05-07", # and this
"245-762a"), # change to this
Serial = factor(Serial))Check that it worked:
temps_doublecheck <- temps %>%
# SAVE the obs for all Serials but 229-070 & 245-762
dplyr::filter((Serial == "229-072" & date_only >= "2021-06-07") |
(Serial == "229-060" & date_only >= "2021-05-18") |
(Serial == "229-074" & date_only >= "2021-04-28") |
(Serial == "229-073" & date_only >= "2021-05-07") |
(Serial == "245-751" & date_only >= "2021-05-07"))
nrow(temps_doublecheck) + nrow(temps2) == nrow(temps)## [1] TRUEtemp_check2 <- temps2 %>%
left_join(check_gone, by = 'Serial') %>%
dplyr::filter(date_only > date_def_gone) %>%
arrange(Serial) %>%
dplyr::select(Serial, date_only, date_def_gone)
temp_check2## [1] Serial date_only date_def_gone
## <0 rows> (or 0-length row.names)## [1] 229-059 229-060 229-063 229-065 229-066 229-069 229-070a 229-070
## [9] 229-072 229-073 229-074 229-077 229-078 229-079 229-080 245-741
## [17] 245-742 245-745 245-748 245-751 245-752 245-753 245-754 245-756
## [25] 245-757 245-759 245-760 245-761 245-762a 245-762 252-881 252-882
## [33] 252-883 252-884 252-885 252-886 252-887
## 37 Levels: 229-059 229-060 229-063 229-065 229-066 229-069 229-070 ... 252-887Subset the temperature data so that we do not use the temperatures from during hydric measurement periods. Also add columns for time interval chunks to make averaging easier later.
# capture weekend dates
cap_dates <- as.Date(c("2021-04-23", "2021-04-24", "2021-04-25", # April
"2021-05-07", "2021-05-08", "2021-05-09"), # May
format = "%Y-%m-%d")
# grouping notes by date
group_dates <- data.frame(date_only = c(seq(as.Date("2021-04-26"),
as.Date("2021-07-13"), by = 1))) %>%
dplyr::filter(date_only %nin% cap_dates) %>%
mutate(date_chunks = factor(c(rep("Apr 26 - May 6", 11),
rep("May 10 - 20", 11),
rep("May 21 - 31", 11),
rep("Jun 1 - 11", 11),
rep("Jun 12 - 22", 11),
rep("Jun 23 - Jul 3", 11),
rep("Jul 4 - 13", 10)),
levels = c("Apr 26 - May 6", "May 10 - 20",
"May 21 - 31", "Jun 1 - 11", "Jun 12 - 22",
"Jun 23 - Jul 3", "Jul 4 - 13")))
# subset dataset and add grouping notes
temp_sub <- temps2 %>%
dplyr::filter(date_only %nin% cap_dates) %>%
left_join(group_dates, by = 'date_only')
summary(temp_sub)## DateTime Serial Freq
## Min. :2021-04-26 00:01:18.00 245-756: 5195 Min. :150.0
## 1st Qu.:2021-05-02 12:14:29.75 245-753: 5183 1st Qu.:150.4
## Median :2021-05-12 02:49:33.50 229-078: 4971 Median :151.1
## Mean :2021-05-15 12:15:37.34 252-886: 4487 Mean :151.0
## 3rd Qu.:2021-05-18 14:59:10.50 229-069: 4457 3rd Qu.:151.5
## Max. :2021-07-13 13:11:42.00 252-881: 4182 Max. :152.0
## (Other):68779
## ppTempC pp date_only hour_of_day
## Min. :-23.99 Min. :1001 Min. :2021-04-26 Min. : 0.00
## 1st Qu.: 32.02 1st Qu.:1561 1st Qu.:2021-05-02 1st Qu.:10.00
## Median : 36.40 Median :1698 Median :2021-05-12 Median :13.00
## Mean : 35.33 Mean :1729 Mean :2021-05-14 Mean :12.61
## 3rd Qu.: 39.14 3rd Qu.:1877 3rd Qu.:2021-05-18 3rd Qu.:16.00
## Max. : 76.93 Max. :2798 Max. :2021-07-13 Max. :23.00
##
## date_chunks
## Apr 26 - May 6:39496
## May 10 - 20 :33790
## May 21 - 31 : 9326
## Jun 1 - 11 : 6163
## Jun 12 - 22 : 3163
## Jun 23 - Jul 3: 2552
## Jul 4 - 13 : 2764I also need to remove outliers, which are pretty standard erroneous points when using temperature-sensitive telemetry with devices external to the animal. I will follow Nicole Gaudenti’s 2021 paper, and exclude “any outliers greater than two standard deviations away from each lizard’s mean Tb”
Look at data distribution:
Get mean and SD for each lizard’s Tb data, then use it to filter out extreme values:
temp_clean <- temp_sub %>%
group_by(Serial) %>%
#summarise(outs = boxplot.stats(ppTempC)$out) # less outliers found this way
mutate(mean_Tb = mean(ppTempC),
sd_Tb = sd(ppTempC)) %>%
mutate(accept_low = mean_Tb - (2*sd_Tb),
accept_high = mean_Tb + (2*sd_Tb)) %>%
dplyr::filter(ppTempC < accept_high & ppTempC > accept_low) %>%
dplyr::select(-mean_Tb, -sd_Tb, -accept_low, -accept_high)
summary(temp_clean)## DateTime Serial Freq
## Min. :2021-04-26 00:01:18.00 245-756: 4984 Min. :150.0
## 1st Qu.:2021-05-02 15:13:54.25 245-753: 4974 1st Qu.:150.4
## Median :2021-05-12 08:52:29.50 229-078: 4686 Median :151.1
## Mean :2021-05-15 18:45:28.49 229-069: 4294 Mean :151.0
## 3rd Qu.:2021-05-22 09:58:02.75 252-886: 4251 3rd Qu.:151.5
## Max. :2021-07-13 13:11:42.00 252-881: 3996 Max. :152.0
## (Other):65645
## ppTempC pp date_only hour_of_day
## Min. : 1.67 Min. :1109 Min. :2021-04-26 Min. : 0.00
## 1st Qu.:32.45 1st Qu.:1564 1st Qu.:2021-05-02 1st Qu.:10.00
## Median :36.48 Median :1694 Median :2021-05-12 Median :13.00
## Mean :35.55 Mean :1721 Mean :2021-05-15 Mean :12.75
## 3rd Qu.:39.09 3rd Qu.:1862 3rd Qu.:2021-05-22 3rd Qu.:16.00
## Max. :53.17 Max. :2798 Max. :2021-07-13 Max. :23.00
##
## date_chunks
## Apr 26 - May 6:36977
## May 10 - 20 :32519
## May 21 - 31 : 9080
## Jun 1 - 11 : 6005
## Jun 12 - 22 : 3067
## Jun 23 - Jul 3: 2494
## Jul 4 - 13 : 2688What proportion of points did we remove?
## [1] 0.045489134.5%
Plot again:
Individual IQRs definitely look reasonable.
Also, save the link between radio frequency and transmitter serial number to help with microhabitat data wrangling.
Microhabitat Use
This data was collected by radio-tracking lizards, finding them, then recording what type of microhabitat they were in. We have to read-in a few files, because at first the formatting was not consistent, then put the dataframes all together.
variables in this dataframe:
- Serial = the radio-collar serial number for that observation (to be used to link to specific individuals)
- DateTime = date and time the temperature was logged
- microhabitat (was Micro_1) = which microhabitat lizards were observed in
- Breed_Col = breeding coloration of lizards (we will not use)
- bpmTempC = the temperature in Celsius from the lizards’ collar, based on the handheld receivers used to locate them, which record the pulse interval in beats per minute rather than the pulse period of the tower receiver for other Tb data
- dwsTempOut = the temperature recorded by a weather station at the study site; these values were interpolated to the specific time of the observation since the weather station was only recording data every two hours ModTempC = the ground surface temperature estimated by a temperature model (Ian Axsom) for the location and time that each observation occurred; this is an estimate of the operative temperature available to these lizards
This is the nice, clean data that was collected early May onwards:
microhab_clean <- read.csv("./data/telemetry/behavior_observations.csv") %>%
# force NAs on, then remove, observations for the "uncollared" lizard
mutate(Serial_no = as.numeric(substr(Serial, 1, 5))) %>%
dplyr::filter(complete.cases(Serial_no)) %>%
# properly format data classes
mutate(date_only = as.Date(substr(DateTime, 1, 10), format = "%Y-%m-%d"),
hour_of_day = as.numeric(substr(DateTime, 12, 13)),
DateTime = as.POSIXct(DateTime, format = "%Y-%m-%d %H:%M:%S"),
Serial = factor(paste(paste(substr(as.character(Serial), 1, 3),
substr(as.character(Serial), 4, 7), sep = "-"))),
microhabitat = factor(Micro_1,
labels = c("Burrow", # keep
"Open", # change from Burrow slide
"Open", # change from Burrow-Partial
"Burrow", # change from Burrow/ephedra
"Burrow", # change from Burrow/shrub
"", ### change from Dead/Missing ###
"Open", # change from Grass trims
"", # change from IDK
"", # change from No Visual
"Open", # keep
"Shade - Full", # keep
"Shade - Partial", # keep
"Shade - Partial", # remove notes
"Shade - Partial", # remove notes
"Shade - Partial", # remove notes
"Shade - Partial", # remove notes
"Shade - Partial", # remove notes
"Shade - Partial", # remove notes
"Shade - Partial" # remove notes
))
) %>%
# for some reason, it only works to fix the order separately
mutate(microhabitat = factor(microhabitat,
levels = c("Open",
"Shade - Partial",
"Shade - Full",
"Burrow"
#"No Visual", "Dead/Missing"
))) %>%
dplyr::select(Serial, DateTime, date_only, hour_of_day, microhabitat,
bpmTempC, dwsTempOut, ModTempC) %>%
# add date interval chunk grouping notes
dplyr::filter(date_only %nin% cap_dates)This is the messy data collected in late April and very early May. Unfortunately, the following radio-frequencies were unable to be linked to a radio collar serial number: 150.14, 150.19, 150.21, and 150.28. But, that only applied to 11 observations.
microhab_early0 <- read.csv("./data/telemetry/behavior_observations_early_0.csv") %>%
dplyr::select(date_only = date, time,
microhabitat = micro,
Serial = serial,
Freq = freq, freq_notes = X)
# this one must be split by whether freq or serial number of radio collars recorded
microhab_early0_freq_only <- microhab_early0 %>%
dplyr::filter(!complete.cases(Serial)) %>%
dplyr::select(-Serial)
microhab_early0_serial <- microhab_early0 %>%
dplyr::filter(complete.cases(Serial)) %>%
mutate(Serial = factor(paste(substr(Serial, 1, 3),
substr(Serial, 4, 7), sep = "-")))
# load other csv, and add the dfs above to it
microhab_early_all <- read.csv("./data/telemetry/behavior_observations_early_1.csv") %>%
dplyr::select(date_only = date, time,
microhabitat = micro,
Freq = freq, freq_notes) %>%
rbind(microhab_early0_freq_only) %>%
mutate(Freq = round(Freq, 2)) %>%
left_join(freq_serial, by = 'Freq') %>%
dplyr::select(-n) %>%
rbind(microhab_early0_serial) %>%
mutate(DateTime = as.POSIXct(paste(date_only, time, sep = " "),
format = "%m/%d/%y %H:%M"),
hour_of_day = as.numeric(substr(as.character(DateTime), 12, 13)),
date_only = as.Date(date_only, format = "%m/%d/%y"),
microhabitat = factor(str_to_title(microhabitat))) %>%
dplyr::filter(complete.cases(Serial, microhabitat)) %>%
# some of these were "focal observations" taken every few minutes
# to be consistent with other obs, take first point of each focal
group_by(Serial, date_only) %>%
#summarise(min(date), max(date), min(time), max(time)) # check that data is one chunk
# might remove a few usable data, but OK
slice(1) %>%
# remove these to make it easy to bind to clean one
dplyr::select(-time, -Freq, -freq_notes) %>%
# add these to make it easy to bind to clean one
mutate(bpmTempC = NA, dwsTempOut = NA, ModTempC = NA)
summary(microhab_early_all) ## date_only microhabitat Serial
## Min. :2021-04-27 Burrow :28 245-752: 5
## 1st Qu.:2021-04-28 Open :79 245-760: 5
## Median :2021-04-29 Shade - Full : 4 252-884: 5
## Mean :2021-04-29 Shade - Partial: 7 229-060: 4
## 3rd Qu.:2021-05-01 229-063: 4
## Max. :2021-05-02 229-069: 4
## (Other):91
## DateTime hour_of_day bpmTempC dwsTempOut
## Min. :2021-05-01 10:03:00.00 Min. : 2.00 Mode:logical Mode:logical
## 1st Qu.:2021-05-01 13:35:30.00 1st Qu.:10.50 NA's:118 NA's:118
## Median :2021-05-01 16:00:00.00 Median :13.00
## Mean :2021-05-01 22:00:40.35 Mean :11.93
## 3rd Qu.:2021-05-02 11:00:30.00 3rd Qu.:15.00
## Max. :2021-05-02 15:46:00.00 Max. :17.00
## NA's :63 NA's :63
## ModTempC
## Mode:logical
## NA's:118
##
##
##
##
## Now, put together the original and messy-now-clean dataframes, and add some formatting:
microhab_pre <- microhab_clean %>%
rbind(microhab_early_all) %>%
dplyr::filter(complete.cases(microhabitat)) %>%
left_join(group_dates, by = 'date_only') %>%
left_join(liz_IDs_only, by = c('Serial' = 'radio_collar_serial')) %>%
mutate(above_below = factor(case_when(microhabitat == "Burrow" ~ "Belowground",
microhabitat != "Burrow" ~ "Aboveground"),
levels = c("Aboveground", "Belowground")))
summary(microhab_pre)## Serial DateTime date_only
## 245-744 : 100 Min. :2021-05-01 10:03:00.00 Min. :2021-04-27
## 245-741 : 96 1st Qu.:2021-05-31 13:30:00.00 1st Qu.:2021-05-29
## 229-078 : 93 Median :2021-06-22 13:17:00.00 Median :2021-06-21
## 245-762a: 92 Mean :2021-06-15 04:09:31.63 Mean :2021-06-12
## 229-079 : 91 3rd Qu.:2021-06-30 14:17:00.00 3rd Qu.:2021-06-30
## 245-745 : 91 Max. :2021-07-13 13:27:00.00 Max. :2021-07-13
## (Other) :1213 NA's :63
## hour_of_day microhabitat bpmTempC dwsTempOut
## Min. : 2.00 Open :544 Min. : 8.028 Min. :10.11
## 1st Qu.:10.00 Shade - Partial:147 1st Qu.:32.828 1st Qu.:26.61
## Median :12.00 Shade - Full :214 Median :36.267 Median :31.28
## Mean :11.99 Burrow :871 Mean :35.633 Mean :30.25
## 3rd Qu.:14.00 3rd Qu.:39.792 3rd Qu.:34.60
## Max. :18.00 Max. :49.823 Max. :40.58
## NA's :63 NA's :150 NA's :173
## ModTempC date_chunks individual_ID sex_M_F
## Min. :21.24 Apr 26 - May 6:166 M-09 : 100 Female: 539
## 1st Qu.:39.50 May 10 - 20 :158 F-10 : 96 Male :1237
## Median :44.48 May 21 - 31 :181 M-14 : 93
## Mean :43.77 Jun 1 - 11 :181 M-03a : 92
## 3rd Qu.:48.99 Jun 12 - 22 :253 M-08 : 91
## Max. :57.86 Jun 23 - Jul 3:579 M-11 : 91
## NA's :118 Jul 4 - 13 :258 (Other):1213
## above_below
## Aboveground:905
## Belowground:871
##
##
##
##
## Check whether we have observations on lizards that were not alive or were not actually wearing their radio-collars because they dropped them.
microhab_check <- microhab_pre %>%
left_join(check_gone, by = 'Serial') %>%
dplyr::filter(date_only > date_def_gone) %>%
arrange(Serial) %>%
dplyr::select(individual_ID, Serial, date_only, date_def_gone)
microhab_check## individual_ID Serial date_only date_def_gone
## 1 F-17 229-060 2021-05-19 2021-05-18
## 2 F-08 229-070 2021-05-11 2021-05-07
## 3 F-08 229-070 2021-05-23 2021-05-07
## 4 F-08 229-070 2021-06-07 2021-05-07
## 5 M-03 245-762 2021-06-10 2021-05-07Check what the name switch should be for F-08 and M-03:
## [1] 245-761 252-886 245-762 252-887 245-753 245-752 245-754 229-079
## [9] 245-744 229-069 252-885 245-759 245-751 245-757 229-059 229-063
## [17] 245-747 229-070 229-065 245-741 245-748 252-882 252-881 229-078
## [25] 245-750 229-080 245-760 245-746 245-745 245-742 252-884 252-883
## [33] 245-756 229-073 229-066 229-072 229-060 229-074 229-077 245-762a
## [41] 229-070a
## 41 Levels: 229-059 229-060 229-063 229-065 229-066 229-069 229-070 ... 252-887## [1] M-01 M-02 M-03 M-04 M-05 M-06 M-07 M-08 M-09 M-10 F-01 F-02
## [13] F-03 F-04 F-05 F-06 F-07 F-08 F-09 F-10 M-11 M-12 M-13 M-14
## [25] M-15 M-16 M-17 M-18 M-19 M-20 F-11 F-12 F-13 F-14 F-15 F-16
## [37] F-17 F-18 F-19 M-03a F-08a
## 79 Levels: F-01 F-02 F-03 F-04 F-05 F-06 F-07 F-08 F-08a F-09 F-10 ... W-039Then fix those things:
microhab <- microhab_pre %>%
# delete the one obs for F-17 that was outside the valid date range
dplyr::filter(!(individual_ID == "F-17" & date_only == "2021-05-19")) %>%
# change the serial and ID of F-08 to be "-a" bc collar was switched
# for F-08a
mutate(individual_ID = replace(individual_ID, # what I want to change
individual_ID == "F-08" & # if this
date_only > "2021-05-07", # and this
"F-08a"), # change to this
Serial = replace(Serial, # what I want to change
Serial == "229-070" & # if this
date_only > "2021-05-07", # and this
"229-070a"), # change to this
# for M-03a
individual_ID = replace(individual_ID, # what I want to change
individual_ID == "M-03" & # if this
date_only > "2021-05-07", # and this
"M-03a"), # change to this
Serial = replace(Serial, # what I want to change
Serial == "245-762" & # if this
date_only > "2021-05-07", # and this
"245-762a") # change to this
)Check that worked:
## [1] TRUEmicrohab_check2 <- microhab %>%
left_join(check_gone, by = 'Serial') %>%
dplyr::filter(date_only > date_def_gone) %>%
arrange(Serial) %>%
dplyr::select(individual_ID, Serial, date_only, date_def_gone)
microhab_check2## [1] individual_ID Serial date_only date_def_gone
## <0 rows> (or 0-length row.names)Export full microhab df
methods stats: How many observations and days was each lizard tracked for?
micro_obs_ns <- microhab %>%
group_by(Serial) %>%
summarise(n = n(),
min_date = min(date_only),
max_date = max(date_only),
date_range = max_date - min_date) %>%
mutate(per_day = n/as.numeric(date_range),
per_wk = per_day*7)
micro_obs_ns %>%
dplyr::filter(Serial != "252-882") %>%
summarise(mean(n), sd(n),
mean(date_range), sd(date_range),
mean(per_wk), sd(per_wk))## # A tibble: 1 × 6
## `mean(n)` `sd(n)` `mean(date_range)` `sd(date_range)` `mean(per_wk)`
## <dbl> <dbl> <drtn> <dbl> <dbl>
## 1 46.7 38.5 47.13158 days 28.8 6.09
## # ℹ 1 more variable: `sd(per_wk)` <dbl>Thermal Metrics
For each lizard, we want to calculate some assessment of how hot they let themselves be, and their thermoregulatory accuracy.
Daily Maximums
daily_max <- temp_clean %>%
# quick look to use
#dplyr::filter(date_chunks %in% c("Apr 26 - May 6", "May 10 - 20")) %>%
# first do daily maximums
group_by(date_only, date_chunks, Serial) %>%
mutate(Tb_percentile_50 = quantile(ppTempC, 0.5),
Tb_percentile_60 = quantile(ppTempC, 0.6),
Tb_percentile_70 = quantile(ppTempC, 0.7),
Tb_percentile_80 = quantile(ppTempC, 0.8),
Tb_percentile_90 = quantile(ppTempC, 0.9)) %>%
ungroup() %>%
# select only the data between the 80-90th percentiles of Tb for each lizard
dplyr::filter(ppTempC > Tb_percentile_80 & ppTempC < Tb_percentile_90) %>%
# calculate the mean Tb for that subset for each lizard per day
group_by(date_only, date_chunks, Serial) %>%
summarise(daily_Tb_percentile_80_90 = mean(ppTempC)) %>%
# then do the mean Tb between 80-90%ile for each lizard per time chunk (~11 days)
group_by(date_chunks, Serial) %>%
summarise(chunk_avg_Tb_percentile_80_90 = mean(daily_Tb_percentile_80_90))
summary(daily_max)## date_chunks Serial chunk_avg_Tb_percentile_80_90
## Apr 26 - May 6:35 229-065: 7 Min. :18.89
## May 10 - 20 :32 229-069: 7 1st Qu.:37.51
## May 21 - 31 :28 229-078: 7 Median :39.38
## Jun 1 - 11 :23 245-742: 7 Mean :38.37
## Jun 12 - 22 :18 245-748: 7 3rd Qu.:40.79
## Jun 23 - Jul 3:15 245-753: 7 Max. :46.85
## Jul 4 - 13 :12 (Other):121Db
From Kat Ivey’s 2020 paper, Tset range is: 32.3±1.2 - 37.5±1.1 deg C.
From Nicole Gaudenti’s 2021 paper, Tset range: 33.2-37.9 deg C.
I will use the widest margin created by both: 32.3 (Ivey) - 37.9 (Gaudenti)
Db <- temp_clean %>%
# we only want daytime Tb for this
# currently using Kat's limits
dplyr::filter(hour_of_day <=19 & hour_of_day >=7) %>%
mutate(Tpref_low = 32.3,
Tpref_high = 37.9) %>%
# calculate Db for every Tb measurement
mutate(Db = case_when(#ppTempC < Tpref_low ~ ppTempC - Tpref_low, # for pos/neg version
ppTempC < Tpref_low ~ Tpref_low - ppTempC, # for abs value version
ppTempC > Tpref_high ~ ppTempC - Tpref_high,
ppTempC < Tpref_high & ppTempC > Tpref_low ~ 0)) %>%
# then avg Db by lizard & date chunk
group_by(date_chunks, Serial) %>%
summarise(chunk_mean_Db_deg_C_diff = mean(Db))
summary(Db)## date_chunks Serial chunk_mean_Db_deg_C_diff
## Apr 26 - May 6:35 229-065: 7 Min. : 0.000
## May 10 - 20 :32 229-069: 7 1st Qu.: 1.033
## May 21 - 31 :29 229-078: 7 Median : 1.532
## Jun 1 - 11 :25 229-079: 7 Mean : 2.453
## Jun 12 - 22 :20 245-742: 7 3rd Qu.: 2.303
## Jun 23 - Jul 3:18 245-745: 7 Max. :16.429
## Jul 4 - 13 :16 (Other):133Microhabitat Use
Proportions of Observations
To calculate microhabitat use, I literally just calculate the proportion of observations for a given individual/time range that were in each microhabitat.
So, for each time chunk, for each lizard (radio-collar serial), I need to get the following values:
- number of total observations
- number of observations in each microhabitat
- proportion of observations in each microhabitat
- aboveground temps based on Ian’s model ?
We do not have a lot of obs for each lizard for every time chunk, so how many observations should there be to make the proportions valid/reasonable? One potential way to do this could be to at least have 1-2 obs for every microhabitat type we are trying to observe (so, minimum points per date chunk either 4 or 8, once I get rid of the irrelevant ones).
# where to cut off n obs?
microhab %>%
# get only 1 per lizard
group_by(date_chunks, Serial) %>%
summarise(n_obs = n()) %>%
ungroup() %>%
# How many lizards have each (n) obs?
group_by(n_obs) %>%
summarise(n = n())## # A tibble: 26 × 2
## n_obs n
## <int> <int>
## 1 1 5
## 2 2 12
## 3 3 6
## 4 4 18
## 5 5 27
## 6 6 18
## 7 7 13
## 8 8 4
## 9 9 7
## 10 10 9
## # ℹ 16 more rowsIt seems like a good cut-off is a minimum of 4 observations per lizard, per time chunk. Otherwise, we would lose a LOT of data.
# micro use BY individual lizard
MH_use_by_indiv <- microhab %>%
group_by(Serial, date_chunks, microhabitat) %>%
summarise(n_obs_micro = n()) %>%
group_by(Serial, date_chunks) %>%
mutate(n_obs_total = sum(n_obs_micro)) %>%
ungroup() %>%
# expect that this takes out 42 obs/rows
dplyr::filter(n_obs_total >= 4) %>%
mutate(MH_use_proportion = n_obs_micro/n_obs_total) #%>%
# check that all proportions for a given lizard and time chunk = 1
# group_by(Serial, date_chunks) %>%
# summarise(prop_total = sum(MH_use_proportion))
head(MH_use_by_indiv)## # A tibble: 6 × 6
## Serial date_chunks microhabitat n_obs_micro n_obs_total MH_use_proportion
## <fct> <fct> <fct> <int> <int> <dbl>
## 1 229-059 Apr 26 - May 6 Open 3 4 0.75
## 2 229-059 Apr 26 - May 6 Burrow 1 4 0.25
## 3 229-060 Apr 26 - May 6 Open 4 6 0.667
## 4 229-060 Apr 26 - May 6 Burrow 2 6 0.333
## 5 229-060 May 10 - 20 Open 2 4 0.5
## 6 229-060 May 10 - 20 Shade - Part… 1 4 0.25## Serial date_chunks microhabitat n_obs_micro
## 245-741 : 23 Apr 26 - May 6:62 Open :122 Min. : 1.000
## 229-079 : 22 May 10 - 20 :67 Shade - Partial: 75 1st Qu.: 1.000
## 245-745 : 22 May 21 - 31 :70 Shade - Full : 67 Median : 3.000
## 245-744 : 21 Jun 1 - 11 :59 Burrow :118 Mean : 4.524
## 245-762a: 19 Jun 12 - 22 :50 3rd Qu.: 5.000
## 245-761 : 18 Jun 23 - Jul 3:45 Max. :32.000
## (Other) :257 Jul 4 - 13 :29
## n_obs_total MH_use_proportion
## Min. : 4.00 Min. :0.0303
## 1st Qu.: 5.00 1st Qu.:0.2000
## Median : 8.00 Median :0.3077
## Mean :11.56 Mean :0.4005
## 3rd Qu.:14.00 3rd Qu.:0.6000
## Max. :37.00 Max. :1.0000
## # micro use for all lizards in the study, aggregated
MH_use_across_indivs <- microhab %>%
group_by(date_chunks, microhabitat) %>%
summarise(n_obs_micro = n(),
mean_mod_temp = mean(ModTempC, na.rm = T)) %>%
group_by(date_chunks) %>%
mutate(n_obs_total = sum(n_obs_micro),
mean_mod_temp2 = mean(mean_mod_temp, na.rm = T)) %>%
ungroup() %>%
mutate(MH_use_proportion = n_obs_micro/n_obs_total)
head(MH_use_across_indivs)## # A tibble: 6 × 7
## date_chunks microhabitat n_obs_micro mean_mod_temp n_obs_total mean_mod_temp2
## <fct> <fct> <int> <dbl> <int> <dbl>
## 1 Apr 26 - Ma… Open 110 37.3 166 40.5
## 2 Apr 26 - Ma… Shade - Par… 13 42.7 166 40.5
## 3 Apr 26 - Ma… Shade - Full 8 41.4 166 40.5
## 4 Apr 26 - Ma… Burrow 35 40.5 166 40.5
## 5 May 10 - 20 Open 82 42.3 157 44.3
## 6 May 10 - 20 Shade - Par… 23 44.3 157 44.3
## # ℹ 1 more variable: MH_use_proportion <dbl># micro use for all lizards, BY SEX
MH_use_across_indivs_sex <- microhab %>%
group_by(date_chunks, microhabitat, sex_M_F) %>%
summarise(n_obs_micro = n(),
mean_mod_temp = mean(ModTempC, na.rm = T)) %>%
group_by(date_chunks, sex_M_F) %>%
mutate(n_obs_total = sum(n_obs_micro),
mean_mod_temp2 = mean(mean_mod_temp, na.rm = T)) %>%
ungroup() %>%
mutate(MH_use_proportion = n_obs_micro/n_obs_total)
summary(MH_use_across_indivs_sex)## date_chunks microhabitat sex_M_F n_obs_micro
## Apr 26 - May 6:8 Open :14 Female:28 Min. : 1.00
## May 10 - 20 :8 Shade - Partial:14 Male :28 1st Qu.: 9.00
## May 21 - 31 :8 Shade - Full :14 Median : 17.00
## Jun 1 - 11 :8 Burrow :14 Mean : 31.70
## Jun 12 - 22 :8 3rd Qu.: 31.75
## Jun 23 - Jul 3:8 Max. :269.00
## Jul 4 - 13 :8
## mean_mod_temp n_obs_total mean_mod_temp2 MH_use_proportion
## Min. :33.22 Min. : 52.0 Min. :39.76 Min. :0.01724
## 1st Qu.:40.21 1st Qu.: 69.0 1st Qu.:40.56 1st Qu.:0.08974
## Median :43.50 Median : 86.5 Median :43.75 Median :0.15281
## Mean :43.08 Mean :126.8 Mean :43.08 Mean :0.25000
## 3rd Qu.:45.85 3rd Qu.:141.0 3rd Qu.:44.92 3rd Qu.:0.37911
## Max. :49.76 Max. :438.0 Max. :45.23 Max. :0.77000
## Proportion of time aboveground
A version of microhabitat use that’s basically burrow/not.
prop_above_by_indiv <- MH_use_by_indiv %>%
dplyr::filter(microhabitat == "Burrow") %>%
mutate(prop_above = 1 - MH_use_proportion)
prop_above_across_indiv <- MH_use_across_indivs %>%
dplyr::filter(microhabitat == "Burrow") %>%
mutate(prop_above = 1 - MH_use_proportion)
prop_above_across_indiv_by_sex <- MH_use_across_indivs_sex %>%
dplyr::filter(microhabitat == "Burrow") %>%
mutate(prop_above = 1 - MH_use_proportion)Onset of Estivation
This didn’t make it into the paper, because the estivation onset dates we got were somewhat questionable, and we have no way of knowing whether it is due to our tracking methods and the dropping/loss of collars throughout the season.
# first, get some key dates for each lizard
last_date <- microhab %>%
# for each lizard
group_by(Serial) %>%
# what was the latest date we tracked them?
summarise(date_last_obs = max(date_only),
# and how many observations did they have total?
n_obs_total = n())
last_date_above <- microhab %>%
# first, subset aboveground obs only
dplyr::filter(microhabitat != "Burrow") %>%
# for each lizard
group_by(Serial) %>%
# what was the latest date we observed each lizard ABOVEground?
summarise(est_onset = max(date_only))
estivating <- microhab %>%
left_join(last_date_above, by = 'Serial') %>%
dplyr::filter(date_only > est_onset) %>%
group_by(Serial, est_onset) %>%
# how many times did we observe each lizard while they were estivating?
summarise(n_obs_est = n()) %>%
left_join(last_date, by = 'Serial') %>%
mutate(est_obs_days = as.numeric(date_last_obs - est_onset),
prop_obs_est = n_obs_est/n_obs_total) %>%
# make sure we tracked them after "last" aboveground observation
dplyr::filter(est_onset < date_last_obs) %>%
# and we still had a signal on them till the end
dplyr::filter(date_last_obs == "2021-07-13")
summary(estivating)## Serial est_onset n_obs_est date_last_obs
## 229-063 :1 Min. :2021-05-11 Min. : 7.00 Min. :2021-07-13
## 229-069 :1 1st Qu.:2021-06-11 1st Qu.:30.25 1st Qu.:2021-07-13
## 229-070a:1 Median :2021-06-20 Median :48.00 Median :2021-07-13
## 229-078 :1 Mean :2021-06-18 Mean :42.29 Mean :2021-07-13
## 229-079 :1 3rd Qu.:2021-06-27 3rd Qu.:58.75 3rd Qu.:2021-07-13
## 245-742 :1 Max. :2021-07-08 Max. :65.00 Max. :2021-07-13
## (Other) :8
## n_obs_total est_obs_days prop_obs_est
## Min. :70.00 Min. : 5.00 Min. :0.07692
## 1st Qu.:83.50 1st Qu.:15.75 1st Qu.:0.38155
## Median :86.50 Median :22.50 Median :0.55868
## Mean :85.57 Mean :24.43 Mean :0.50250
## 3rd Qu.:90.75 3rd Qu.:31.25 3rd Qu.:0.65830
## Max. :93.00 Max. :63.00 Max. :0.92857
## estivating %>%
mutate(est = "Yes") %>%
group_by(est) %>%
summarise(mean(n_obs_est),
sd(n_obs_est),
min(n_obs_est),
max(n_obs_est))## # A tibble: 1 × 5
## est `mean(n_obs_est)` `sd(n_obs_est)` `min(n_obs_est)` `max(n_obs_est)`
## <chr> <dbl> <dbl> <int> <int>
## 1 Yes 42.3 19.8 7 65Max Temps ~ Hydric
Now, look for relationships between hydric physiology and thermoregulation.
Prep Data
max_temps_hydric <- liz_dat_apr_may_only %>%
left_join(daily_max, by = c('date_chunks',
'radio_collar_serial' = 'Serial')) %>%
# remove NAs
dplyr::filter(complete.cases(chunk_avg_Tb_percentile_80_90)) %>%
# remove the one outlier wayyyyy below the other values
dplyr::filter(chunk_avg_Tb_percentile_80_90>30)
summary(max_temps_hydric)## capture_date_time capture_date radio_collar_serial
## Min. :2021-04-23 10:39:00 Min. :2021-04-23 229-059: 2
## 1st Qu.:2021-04-23 14:06:30 1st Qu.:2021-04-23 229-060: 2
## Median :2021-04-24 10:49:00 Median :2021-04-24 229-063: 2
## Mean :2021-04-28 10:39:18 Mean :2021-04-28 229-069: 2
## 3rd Qu.:2021-05-07 11:27:00 3rd Qu.:2021-05-07 229-072: 2
## Max. :2021-05-08 14:19:00 Max. :2021-05-08 229-077: 2
## NA's :4 (Other):42
## PIT_tag_ID individual_ID sex_M_F mass_g SVL_mm
## Length:54 F-01 : 2 Female:27 Min. :26.00 Min. : 85.0
## Class :character F-05 : 2 Male :27 1st Qu.:33.17 1st Qu.: 95.0
## Mode :character F-06 : 2 Median :37.50 Median : 99.0
## F-11 : 2 Mean :38.93 Mean :100.4
## F-12 : 2 3rd Qu.:43.75 3rd Qu.:107.8
## F-13 : 2 Max. :53.60 Max. :122.0
## (Other):42
## tmt tmt_mass_change_g capture_to_msmt hematocrit_percent
## Water Tmt:25 Min. :-3.0000 Min. : 11.00 Min. :23.00
## Sham Tmt :29 1st Qu.:-1.0000 1st Qu.: 59.75 1st Qu.:30.00
## No Tmt : 0 Median : 0.0000 Median : 152.00 Median :36.00
## Mean :-0.3111 Mean : 410.56 Mean :34.81
## 3rd Qu.: 0.0000 3rd Qu.: 265.00 3rd Qu.:39.00
## Max. : 3.0000 Max. :1665.00 Max. :58.00
## NA's :4 NA's :1
## cloacal_temp_C CEWL_g_m2h msmt_temp_C msmt_RH_percent
## Min. :24.0 Min. : 0.650 Min. :19.07 Min. :12.20
## 1st Qu.:32.0 1st Qu.: 7.077 1st Qu.:28.71 1st Qu.:14.53
## Median :33.0 Median : 9.262 Median :30.62 Median :17.87
## Mean :32.5 Mean : 9.481 Mean :29.49 Mean :19.30
## 3rd Qu.:34.0 3rd Qu.:11.377 3rd Qu.:31.86 3rd Qu.:22.38
## Max. :36.0 Max. :21.673 Max. :33.55 Max. :35.83
##
## e_s_kPa msmt_VPD_kPa osmolality_mmol_kg_mean month week
## Min. :2.205 Min. :1.415 Min. :316.0 April:34 1 :34
## 1st Qu.:3.938 1st Qu.:3.064 1st Qu.:347.7 May :20 3 :20
## Median :4.396 Median :3.694 Median :365.0 July : 0 12: 0
## Mean :4.188 Mean :3.421 Mean :367.3
## 3rd Qu.:4.718 3rd Qu.:4.035 3rd Qu.:381.3
## Max. :5.188 Max. :4.319 Max. :436.5
## NA's :1
## week_num ultrasound_date gravid_Y_N number_eggs
## Min. :1.000 Min. :2021-04-24 Gravid :14 Min. :3.000
## 1st Qu.:1.000 1st Qu.:2021-04-24 Not Gravid:11 1st Qu.:3.250
## Median :1.000 Median :2021-04-24 NA's :29 Median :4.000
## Mean :1.741 Mean :2021-04-29 Mean :3.714
## 3rd Qu.:3.000 3rd Qu.:2021-05-08 3rd Qu.:4.000
## Max. :3.000 Max. :2021-05-08 Max. :4.000
## NA's :29 NA's :40
## largest_egg_diameter_cm largest_egg_circumference_cm stage
## Min. :0.6400 Min. :1.820 small round: 1
## 1st Qu.:0.8775 1st Qu.:2.502 large round: 8
## Median :0.9950 Median :2.885 soft : 2
## Mean :1.1093 Mean :3.140 soft oblong: 0
## 3rd Qu.:1.3750 3rd Qu.:3.970 firm oblong: 2
## Max. :1.8000 Max. :4.740 hard oblong: 1
## NA's :40 NA's :40 NA's :40
## dev_point SMI date_chunks
## Min. :1.000 Min. :27.04 Apr 26 - May 6:34
## 1st Qu.:2.000 1st Qu.:33.00 May 10 - 20 :20
## Median :2.000 Median :34.73 May 21 - 31 : 0
## Mean :2.714 Mean :35.52 Jun 1 - 11 : 0
## 3rd Qu.:3.000 3rd Qu.:38.44 Jun 12 - 22 : 0
## Max. :5.000 Max. :45.63 Jun 23 - Jul 3: 0
## NA's :40 Jul 4 - 13 : 0
## chunk_avg_Tb_percentile_80_90
## Min. :35.53
## 1st Qu.:38.35
## Median :39.14
## Mean :39.05
## 3rd Qu.:39.81
## Max. :41.64
## Plots
sub_maxs <- max_temps_hydric %>%
dplyr::select(chunk_avg_Tb_percentile_80_90, CEWL_g_m2h, osmolality_mmol_kg_mean, mass_g, SMI)
PerformanceAnalytics::chart.Correlation(sub_maxs)
ggplot(max_temps_hydric) +
geom_point(aes(x = CEWL_g_m2h,
y = chunk_avg_Tb_percentile_80_90,
color = date_chunks)) +
geom_smooth(aes(x = CEWL_g_m2h,
y = chunk_avg_Tb_percentile_80_90),
se = F,
method = "lm",
formula = y~x)
ggplot(max_temps_hydric) +
geom_point(aes(x = osmolality_mmol_kg_mean,
y = chunk_avg_Tb_percentile_80_90,
color = date_chunks)) +
geom_smooth(aes(x = osmolality_mmol_kg_mean,
y = chunk_avg_Tb_percentile_80_90),
se = F,
method = "lm",
formula = y~x)
Models
# CEWL
max_temp_CEWL_mod <- lmerTest::lmer(data = max_temps_hydric,
chunk_avg_Tb_percentile_80_90 ~ CEWL_g_m2h +
(1|individual_ID))
summary(max_temp_CEWL_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: chunk_avg_Tb_percentile_80_90 ~ CEWL_g_m2h + (1 | individual_ID)
## Data: max_temps_hydric
##
## REML criterion at convergence: 174.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.07592 -0.47459 -0.06368 0.51363 1.69421
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID (Intercept) 0.6737 0.8208
## Residual 0.7941 0.8911
## Number of obs: 54, groups: individual_ID, 36
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 39.13941 0.38430 48.59891 101.845 <2e-16 ***
## CEWL_g_m2h -0.02071 0.03586 39.07870 -0.578 0.567
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## CEWL_g_m2h -0.876
# osmolality
max_temp_osml_mod <- lmerTest::lmer(data = max_temps_hydric,
chunk_avg_Tb_percentile_80_90 ~ osmolality_mmol_kg_mean +
(1|individual_ID))
summary(max_temp_osml_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: chunk_avg_Tb_percentile_80_90 ~ osmolality_mmol_kg_mean + (1 |
## individual_ID)
## Data: max_temps_hydric
##
## REML criterion at convergence: 172.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0620 -0.5175 -0.0073 0.5198 1.6943
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID (Intercept) 0.7480 0.8649
## Residual 0.6931 0.8325
## Number of obs: 53, groups: individual_ID, 35
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 35.262208 2.058746 45.768535 17.128 <2e-16 ***
## osmolality_mmol_kg_mean 0.010069 0.005586 45.564647 1.803 0.0781 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## osmllty_m__ -0.996
# mass
max_temp_mass_mod <- lmerTest::lmer(data = max_temps_hydric,
chunk_avg_Tb_percentile_80_90 ~ mass_g +
(1|individual_ID))
summary(max_temp_mass_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: chunk_avg_Tb_percentile_80_90 ~ mass_g + (1 | individual_ID)
## Data: max_temps_hydric
##
## REML criterion at convergence: 173.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.16856 -0.44986 -0.05865 0.47082 1.81093
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID (Intercept) 0.6622 0.8138
## Residual 0.7724 0.8789
## Number of obs: 54, groups: individual_ID, 36
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 37.83524 0.89481 40.51821 42.283 <2e-16 ***
## mass_g 0.02870 0.02266 40.04802 1.267 0.213
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## mass_g -0.979
# SMI
max_temp_SMI_mod <- lmerTest::lmer(data = max_temps_hydric,
chunk_avg_Tb_percentile_80_90 ~ SMI +
(1|individual_ID))
summary(max_temp_SMI_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: chunk_avg_Tb_percentile_80_90 ~ SMI + (1 | individual_ID)
## Data: max_temps_hydric
##
## REML criterion at convergence: 169.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.15402 -0.50939 -0.09287 0.40471 1.99940
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID (Intercept) 0.5610 0.7490
## Residual 0.7738 0.8797
## Number of obs: 54, groups: individual_ID, 36
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 36.18079 1.29543 40.90388 27.93 <2e-16 ***
## SMI 0.07826 0.03623 40.50454 2.16 0.0368 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## SMI -0.991
# supplemental hydration tmt
max_temp_tmt_mod <- lmerTest::lmer(data = max_temps_hydric,
chunk_avg_Tb_percentile_80_90 ~ tmt +
(1|individual_ID))
summary(max_temp_tmt_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: chunk_avg_Tb_percentile_80_90 ~ tmt + (1 | individual_ID)
## Data: max_temps_hydric
##
## REML criterion at convergence: 166.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0046 -0.5040 -0.0538 0.5567 1.6822
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID (Intercept) 0.5918 0.7693
## Residual 0.7909 0.8893
## Number of obs: 54, groups: individual_ID, 36
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 38.6460 0.2542 35.9445 152.052 <2e-16 ***
## tmtSham Tmt 0.6050 0.3578 31.6304 1.691 0.101
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## tmtSham Tmt -0.710
Db ~ Hydric
Prep Data
db_hydric <- liz_dat_apr_may_only %>%
left_join(Db, by = c('date_chunks',
'radio_collar_serial' = 'Serial')) %>%
# remove NAs
dplyr::filter(complete.cases(chunk_mean_Db_deg_C_diff)) %>%
# remove outliers
dplyr::filter(chunk_mean_Db_deg_C_diff > -4)
summary(db_hydric)## capture_date_time capture_date radio_collar_serial
## Min. :2021-04-23 10:39:00.00 Min. :2021-04-23 229-059: 2
## 1st Qu.:2021-04-23 14:13:00.00 1st Qu.:2021-04-23 229-060: 2
## Median :2021-04-24 10:49:00.00 Median :2021-04-24 229-063: 2
## Mean :2021-04-28 08:49:02.34 Mean :2021-04-28 229-069: 2
## 3rd Qu.:2021-05-07 11:23:00.00 3rd Qu.:2021-05-07 229-072: 2
## Max. :2021-05-08 14:19:00.00 Max. :2021-05-08 229-077: 2
## NA's :4 (Other):43
## PIT_tag_ID individual_ID sex_M_F mass_g SVL_mm
## Length:55 F-01 : 2 Female:28 Min. :26.00 Min. : 85.0
## Class :character F-05 : 2 Male :27 1st Qu.:33.35 1st Qu.: 95.0
## Mode :character F-06 : 2 Median :38.00 Median : 99.0
## F-11 : 2 Mean :38.91 Mean :100.6
## F-12 : 2 3rd Qu.:43.60 3rd Qu.:108.0
## F-13 : 2 Max. :53.60 Max. :122.0
## (Other):43
## tmt tmt_mass_change_g capture_to_msmt hematocrit_percent
## Water Tmt:26 Min. :-3.0000 Min. : 11.0 Min. :23.00
## Sham Tmt :29 1st Qu.:-1.0000 1st Qu.: 60.5 1st Qu.:30.00
## No Tmt : 0 Median : 0.0000 Median : 150.0 Median :36.00
## Mean :-0.3236 Mean : 403.8 Mean :34.63
## 3rd Qu.: 0.0000 3rd Qu.: 264.0 3rd Qu.:38.75
## Max. : 3.0000 Max. :1665.0 Max. :58.00
## NA's :4 NA's :1
## cloacal_temp_C CEWL_g_m2h msmt_temp_C msmt_RH_percent
## Min. :24.00 Min. : 0.650 Min. :19.07 Min. :12.20
## 1st Qu.:32.00 1st Qu.: 7.377 1st Qu.:28.74 1st Qu.:14.53
## Median :33.00 Median : 9.318 Median :30.68 Median :18.67
## Mean :32.53 Mean : 9.649 Mean :29.54 Mean :19.32
## 3rd Qu.:34.00 3rd Qu.:11.886 3rd Qu.:31.88 3rd Qu.:22.30
## Max. :36.00 Max. :21.673 Max. :33.55 Max. :35.83
##
## e_s_kPa msmt_VPD_kPa osmolality_mmol_kg_mean month week
## Min. :2.205 Min. :1.415 Min. :316.0 April:35 1 :35
## 1st Qu.:3.945 1st Qu.:3.068 1st Qu.:346.0 May :20 3 :20
## Median :4.411 Median :3.703 Median :364.2 July : 0 12: 0
## Mean :4.199 Mean :3.428 Mean :366.9
## 3rd Qu.:4.723 3rd Qu.:4.033 3rd Qu.:381.2
## Max. :5.188 Max. :4.319 Max. :436.5
## NA's :1
## week_num ultrasound_date gravid_Y_N number_eggs
## Min. :1.000 Min. :2021-04-24 Gravid :14 Min. :3.000
## 1st Qu.:1.000 1st Qu.:2021-04-24 Not Gravid:12 1st Qu.:3.250
## Median :1.000 Median :2021-04-24 NA's :29 Median :4.000
## Mean :1.727 Mean :2021-04-29 Mean :3.714
## 3rd Qu.:3.000 3rd Qu.:2021-05-08 3rd Qu.:4.000
## Max. :3.000 Max. :2021-05-08 Max. :4.000
## NA's :29 NA's :41
## largest_egg_diameter_cm largest_egg_circumference_cm stage
## Min. :0.6400 Min. :1.820 small round: 1
## 1st Qu.:0.8775 1st Qu.:2.502 large round: 8
## Median :0.9950 Median :2.885 soft : 2
## Mean :1.1093 Mean :3.140 soft oblong: 0
## 3rd Qu.:1.3750 3rd Qu.:3.970 firm oblong: 2
## Max. :1.8000 Max. :4.740 hard oblong: 1
## NA's :41 NA's :41 NA's :41
## dev_point SMI date_chunks chunk_mean_Db_deg_C_diff
## Min. :1.000 Min. :27.04 Apr 26 - May 6:35 Min. : 0.5738
## 1st Qu.:2.000 1st Qu.:33.00 May 10 - 20 :20 1st Qu.: 1.0201
## Median :2.000 Median :34.71 May 21 - 31 : 0 Median : 1.2748
## Mean :2.714 Mean :35.37 Jun 1 - 11 : 0 Mean : 1.7030
## 3rd Qu.:3.000 3rd Qu.:38.43 Jun 12 - 22 : 0 3rd Qu.: 1.6024
## Max. :5.000 Max. :45.63 Jun 23 - Jul 3: 0 Max. :11.5644
## NA's :41 Jul 4 - 13 : 0Plots
sub_db <- db_hydric %>%
dplyr::select(chunk_mean_Db_deg_C_diff, CEWL_g_m2h, osmolality_mmol_kg_mean, mass_g, SMI)
PerformanceAnalytics::chart.Correlation(sub_db)
ggplot(db_hydric) +
geom_point(aes(x = CEWL_g_m2h,
y = chunk_mean_Db_deg_C_diff,
color = date_chunks)) +
geom_smooth(aes(x = CEWL_g_m2h,
y = chunk_mean_Db_deg_C_diff),
se = F,
method = "lm",
formula = y~x)
ggplot(db_hydric) +
geom_point(aes(x = osmolality_mmol_kg_mean,
y = chunk_mean_Db_deg_C_diff,
color = date_chunks)) +
geom_smooth(aes(x = osmolality_mmol_kg_mean,
y = chunk_mean_Db_deg_C_diff),
se = F,
method = "lm",
formula = y~x)
Models
# CEWL
db_CEWL_mod <- lmerTest::lmer(data = db_hydric,
chunk_mean_Db_deg_C_diff ~ CEWL_g_m2h +
(1|individual_ID))
summary(db_CEWL_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: chunk_mean_Db_deg_C_diff ~ CEWL_g_m2h + (1 | individual_ID)
## Data: db_hydric
##
## REML criterion at convergence: 193.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.62001 -0.24244 -0.01936 0.31381 1.60181
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID (Intercept) 5.71708 2.3910
## Residual 0.08465 0.2909
## Number of obs: 55, groups: individual_ID, 37
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.10721 0.42177 44.35973 4.996 9.6e-06 ***
## CEWL_g_m2h -0.02517 0.01539 17.28891 -1.635 0.12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## CEWL_g_m2h -0.349
# osmolality
db_osml_mod <- lmerTest::lmer(data = db_hydric,
chunk_mean_Db_deg_C_diff ~ osmolality_mmol_kg_mean +
(1|individual_ID))
summary(db_osml_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## chunk_mean_Db_deg_C_diff ~ osmolality_mmol_kg_mean + (1 | individual_ID)
## Data: db_hydric
##
## REML criterion at convergence: 189.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.20314 -0.29959 -0.02697 0.34570 1.15461
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID (Intercept) 5.689 2.3852
## Residual 0.072 0.2683
## Number of obs: 54, groups: individual_ID, 36
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.634993 1.016775 24.263771 -0.625 0.5381
## osmolality_mmol_kg_mean 0.006915 0.002553 17.693153 2.708 0.0146 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## osmllty_m__ -0.920
# mass
db_mass_mod <- lmerTest::lmer(data = db_hydric,
chunk_mean_Db_deg_C_diff ~ mass_g +
(1|individual_ID))
summary(db_mass_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: chunk_mean_Db_deg_C_diff ~ mass_g + (1 | individual_ID)
## Data: db_hydric
##
## REML criterion at convergence: 192.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.57300 -0.13496 -0.02122 0.17189 1.58553
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID (Intercept) 5.73934 2.3957
## Residual 0.08222 0.2867
## Number of obs: 55, groups: individual_ID, 37
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.42166 0.90373 33.57494 0.467 0.6438
## mass_g 0.03762 0.02115 23.77678 1.779 0.0881 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## mass_g -0.899
# SMI
db_SMI_mod <- lmerTest::lmer(data = db_hydric,
chunk_mean_Db_deg_C_diff ~ SMI +
(1|individual_ID))
summary(db_SMI_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: chunk_mean_Db_deg_C_diff ~ SMI + (1 | individual_ID)
## Data: db_hydric
##
## REML criterion at convergence: 194.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5647 -0.1435 -0.0334 0.1880 1.5731
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID (Intercept) 5.63424 2.3737
## Residual 0.09294 0.3049
## Number of obs: 55, groups: individual_ID, 37
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.21859 0.69212 33.91620 1.761 0.0873 .
## SMI 0.01846 0.01624 17.27587 1.137 0.2712
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## SMI -0.824
# supplemental water tmt
db_tmt_mod <- lmerTest::lmer(data = db_hydric,
chunk_mean_Db_deg_C_diff ~ tmt +
(1|individual_ID))
summary(db_tmt_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: chunk_mean_Db_deg_C_diff ~ tmt + (1 | individual_ID)
## Data: db_hydric
##
## REML criterion at convergence: 187.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.45783 -0.18050 -0.03887 0.31549 1.48654
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID (Intercept) 5.55610 2.3571
## Residual 0.09779 0.3127
## Number of obs: 55, groups: individual_ID, 37
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.1925 0.5310 34.8339 4.129 0.000216 ***
## tmtSham Tmt -0.7086 0.7826 34.6996 -0.905 0.371490
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## tmtSham Tmt -0.678Microhabitat ~ Hydric
Prep Data
I actually need to prep the microhabitat data- the metric I will investigate is proportion of time aboveground.
prop_above_hydric <- prop_above_by_indiv %>%
dplyr::filter(date_chunks %in% c("Apr 26 - May 6", "May 10 - 20"))micro_hydric <- liz_dat_apr_may_only %>%
left_join(prop_above_hydric, by = c('date_chunks',
'radio_collar_serial' = 'Serial')) %>%
# remove NAs
dplyr::filter(complete.cases(prop_above))
summary(micro_hydric)## capture_date_time capture_date radio_collar_serial
## Min. :2021-04-23 10:43:00.00 Min. :2021-04-23 229-060: 2
## 1st Qu.:2021-04-23 13:44:00.00 1st Qu.:2021-04-23 252-884: 2
## Median :2021-04-24 10:49:00.00 Median :2021-04-24 252-887: 2
## Mean :2021-04-28 07:13:02.60 Mean :2021-04-28 229-059: 1
## 3rd Qu.:2021-05-07 11:15:30.00 3rd Qu.:2021-05-08 229-066: 1
## Max. :2021-05-08 14:19:00.00 Max. :2021-05-08 229-070: 1
## NA's :2 (Other):16
## PIT_tag_ID individual_ID sex_M_F mass_g SVL_mm
## Length:25 F-11 : 2 Female:15 Min. :26.00 Min. : 85.0
## Class :character F-17 : 2 Male :10 1st Qu.:32.00 1st Qu.: 94.0
## Mode :character M-04 : 2 Median :35.60 Median : 97.0
## F-01 : 1 Mean :36.81 Mean : 98.4
## F-03 : 1 3rd Qu.:38.00 3rd Qu.:104.0
## F-05 : 1 Max. :52.70 Max. :114.0
## (Other):16
## tmt tmt_mass_change_g capture_to_msmt hematocrit_percent
## Water Tmt:11 Min. :-3.000 Min. : 25.0 Min. :23.00
## Sham Tmt :14 1st Qu.:-1.000 1st Qu.: 67.5 1st Qu.:27.75
## No Tmt : 0 Median : 0.000 Median : 154.0 Median :32.50
## Mean :-0.464 Mean : 437.7 Mean :34.46
## 3rd Qu.: 0.000 3rd Qu.: 267.5 3rd Qu.:40.25
## Max. : 1.800 Max. :1662.0 Max. :58.00
## NA's :2 NA's :1
## cloacal_temp_C CEWL_g_m2h msmt_temp_C msmt_RH_percent
## Min. :24.00 Min. : 0.650 Min. :19.07 Min. :13.12
## 1st Qu.:32.00 1st Qu.: 8.457 1st Qu.:28.46 1st Qu.:14.63
## Median :33.00 Median : 9.495 Median :30.68 Median :16.30
## Mean :32.68 Mean :10.365 Mean :29.10 Mean :19.49
## 3rd Qu.:34.00 3rd Qu.:12.213 3rd Qu.:31.80 3rd Qu.:22.15
## Max. :36.00 Max. :21.673 Max. :33.07 Max. :35.83
##
## e_s_kPa msmt_VPD_kPa osmolality_mmol_kg_mean month week
## Min. :2.205 Min. :1.415 Min. :316.0 April:16 1 :16
## 1st Qu.:3.881 1st Qu.:2.990 1st Qu.:340.2 May : 9 3 : 9
## Median :4.411 Median :3.703 Median :357.8 July : 0 12: 0
## Mean :4.109 Mean :3.358 Mean :359.4
## 3rd Qu.:4.701 3rd Qu.:4.007 3rd Qu.:378.6
## Max. :5.049 Max. :4.309 Max. :402.0
## NA's :1
## week_num ultrasound_date gravid_Y_N number_eggs
## Min. :1.00 Min. :2021-04-24 Gravid : 9 Min. :3.000
## 1st Qu.:1.00 1st Qu.:2021-04-24 Not Gravid: 5 1st Qu.:3.000
## Median :1.00 Median :2021-05-01 NA's :11 Median :4.000
## Mean :1.72 Mean :2021-05-01 Mean :3.667
## 3rd Qu.:3.00 3rd Qu.:2021-05-08 3rd Qu.:4.000
## Max. :3.00 Max. :2021-05-08 Max. :4.000
## NA's :11 NA's :16
## largest_egg_diameter_cm largest_egg_circumference_cm stage
## Min. :0.640 Min. :1.82 small round: 1
## 1st Qu.:0.870 1st Qu.:2.45 large round: 6
## Median :0.990 Median :2.88 soft : 1
## Mean :1.043 Mean :2.95 soft oblong: 0
## 3rd Qu.:1.200 3rd Qu.:3.48 firm oblong: 1
## Max. :1.510 Max. :4.07 hard oblong: 0
## NA's :16 NA's :16 NA's :16
## dev_point SMI date_chunks microhabitat
## Min. :1.000 Min. :28.35 Apr 26 - May 6:16 Open : 0
## 1st Qu.:2.000 1st Qu.:32.48 May 10 - 20 : 9 Shade - Partial: 0
## Median :2.000 Median :34.49 May 21 - 31 : 0 Shade - Full : 0
## Mean :2.333 Mean :35.16 Jun 1 - 11 : 0 Burrow :25
## 3rd Qu.:2.000 3rd Qu.:36.44 Jun 12 - 22 : 0
## Max. :5.000 Max. :45.63 Jun 23 - Jul 3: 0
## NA's :16 Jul 4 - 13 : 0
## n_obs_micro n_obs_total MH_use_proportion prop_above
## Min. :1.00 Min. :4.0 Min. :0.1429 Min. :0.2500
## 1st Qu.:1.00 1st Qu.:4.0 1st Qu.:0.2000 1st Qu.:0.6000
## Median :1.00 Median :5.0 Median :0.2500 Median :0.7500
## Mean :1.68 Mean :5.2 Mean :0.3287 Mean :0.6713
## 3rd Qu.:2.00 3rd Qu.:6.0 3rd Qu.:0.4000 3rd Qu.:0.8000
## Max. :4.00 Max. :7.0 Max. :0.7500 Max. :0.8571
## Plots
sub_micro <- micro_hydric %>%
dplyr::select(prop_above, CEWL_g_m2h, osmolality_mmol_kg_mean, mass_g, SMI)
PerformanceAnalytics::chart.Correlation(sub_micro)
ggplot(MH_use_by_indiv) +
geom_bar(aes(x = date_chunks,
y = MH_use_proportion,
fill = microhabitat
),
position = "fill",
stat = "identity") +
facet_wrap(~Serial) +
theme_classic()
ggplot(micro_hydric) +
geom_point(aes(x = CEWL_g_m2h,
y = prop_above,
color = date_chunks)) +
geom_smooth(aes(x = CEWL_g_m2h,
y = prop_above),
se = F,
method = "lm",
formula = y~x)
ggplot(micro_hydric) +
geom_point(aes(x = osmolality_mmol_kg_mean,
y = prop_above,
color = date_chunks)) +
geom_smooth(aes(x = osmolality_mmol_kg_mean,
y = prop_above),
se = F,
method = "lm",
formula = y~x)
Models
# CEWL
micro_CEWL_mod <- lmerTest::lmer(data = micro_hydric,
prop_above ~ CEWL_g_m2h +
(1|individual_ID))
summary(micro_CEWL_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: prop_above ~ CEWL_g_m2h + (1 | individual_ID)
## Data: micro_hydric
##
## REML criterion at convergence: -9.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.99464 -0.18300 0.07242 0.16112 0.99194
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID (Intercept) 0.03197 0.17881
## Residual 0.00143 0.03781
## Number of obs: 25, groups: individual_ID, 22
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.791099 0.059783 7.601336 13.23 1.61e-06 ***
## CEWL_g_m2h -0.012504 0.004465 2.748072 -2.80 0.075 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## CEWL_g_m2h -0.759
# osmolality
micro_osml_mod <- lmerTest::lmer(data = micro_hydric,
prop_above ~ osmolality_mmol_kg_mean +
(1|individual_ID))
summary(micro_osml_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: prop_above ~ osmolality_mmol_kg_mean + (1 | individual_ID)
## Data: micro_hydric
##
## REML criterion at convergence: -3.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.0054 -0.4661 0.1346 0.3529 0.9514
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID (Intercept) 0.026216 0.16191
## Residual 0.005785 0.07606
## Number of obs: 24, groups: individual_ID, 21
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.031922 0.420963 8.191928 2.451 0.0392 *
## osmolality_mmol_kg_mean -0.001043 0.001166 8.075337 -0.895 0.3966
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## osmllty_m__ -0.996
# mass
micro_mass_mod <- lmerTest::lmer(data = micro_hydric,
prop_above ~ mass_g +
(1|individual_ID))
summary(micro_mass_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: prop_above ~ mass_g + (1 | individual_ID)
## Data: micro_hydric
##
## REML criterion at convergence: -6.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.98776 -0.41604 0.06801 0.40571 0.95433
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID (Intercept) 0.025796 0.16061
## Residual 0.006184 0.07864
## Number of obs: 25, groups: individual_ID, 22
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.524904 0.187152 21.190865 2.805 0.0106 *
## mass_g 0.003797 0.004978 21.259937 0.763 0.4540
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## mass_g -0.979
# SMI
micro_SMI_mod <- lmerTest::lmer(data = micro_hydric,
prop_above ~ SMI +
(1|individual_ID))
summary(micro_SMI_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: prop_above ~ SMI + (1 | individual_ID)
## Data: micro_hydric
##
## REML criterion at convergence: -7.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.0319 -0.3658 0.1095 0.3462 0.9900
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID (Intercept) 0.028661 0.16930
## Residual 0.004752 0.06893
## Number of obs: 25, groups: individual_ID, 22
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.472810 0.306820 22.773223 1.541 0.137
## SMI 0.005427 0.008619 22.720145 0.630 0.535
##
## Correlation of Fixed Effects:
## (Intr)
## SMI -0.992
# supplemental hydration tmt
micro_tmt_mod <- lmerTest::lmer(data = micro_hydric,
prop_above ~ tmt +
(1|individual_ID))
summary(micro_tmt_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: prop_above ~ tmt + (1 | individual_ID)
## Data: micro_hydric
##
## REML criterion at convergence: -12.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.0381 -0.2972 0.1325 0.3188 0.9215
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID (Intercept) 0.027290 0.1652
## Residual 0.005315 0.0729
## Number of obs: 25, groups: individual_ID, 22
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.63293 0.05685 19.84699 11.133 5.52e-10 ***
## tmtSham Tmt 0.05782 0.07687 19.72667 0.752 0.461
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## tmtSham Tmt -0.740All Anovas
Just so they’re all in one place, get all the anovas from the models above to run:
CEWL
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## CEWL_g_m2h 0.24422 0.24422 1 42.422 0.3076 0.5821## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## CEWL_g_m2h 0.22482 0.22482 1 17.734 2.6559 0.1208## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## CEWL_g_m2h 0.008691 0.008691 1 2.8375 6.0795 0.09528 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Osmolality
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## osmolality_mmol_kg_mean 2.073 2.073 1 46.788 2.9908 0.09033 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## osmolality_mmol_kg_mean 0.52367 0.52367 1 17.855 7.2737 0.01482 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## osmolality_mmol_kg_mean 0.0030997 0.0030997 1 7.7022 0.5358 0.4858Body Mass
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## mass_g 1.1961 1.1961 1 41.917 1.5485 0.2203## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## mass_g 0.24808 0.24808 1 24.158 3.0171 0.09513 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## mass_g 0.0034519 0.0034519 1 21.557 0.5581 0.4631Body Condition
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## SMI 3.3291 3.3291 1 41.314 4.3023 0.04433 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## SMI 0.11927 0.11927 1 17.676 1.2833 0.2724## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## SMI 0.0015565 0.0015565 1 22.753 0.3276 0.5727Correlation “Matrix”
For each model, get the marginal Rsq for the amount of variance that’s explained by fixed effects only. Then, we will use that to calculate R, and make a figure. I put “matrix” in quotations because it will be a subset of what would be presented for a full correlation matrix.
corrs <- data.frame(hydric_var = c(rep("CEWL", 3), rep("Plasma Osmolality", 3),
#rep("Body Mass", 3),
rep("Body Condition", 3)),
thermal_var = c(rep(c("Max Temp", "Db", "Micro Use"), 3)),
R = c(
# CEWL
cor(max_temps_hydric$CEWL_g_m2h,
max_temps_hydric$chunk_avg_Tb_percentile_80_90),
cor(db_hydric$CEWL_g_m2h, db_hydric$chunk_mean_Db_deg_C_diff),
cor(micro_hydric$CEWL_g_m2h, micro_hydric$prop_above),
# osmolality
cor(max_temps_hydric$osmolality_mmol_kg_mean,
max_temps_hydric$chunk_avg_Tb_percentile_80_90,
use = "complete.obs"),
cor(db_hydric$osmolality_mmol_kg_mean,
db_hydric$chunk_mean_Db_deg_C_diff,
use = "complete.obs"),
cor(micro_hydric$osmolality_mmol_kg_mean, micro_hydric$prop_above,
use = "complete.obs"),
# mass
# cor(max_temps_hydric$mass_g,
# max_temps_hydric$chunk_avg_Tb_percentile_80_90),
# cor(db_hydric$mass_g, db_hydric$chunk_mean_Db_deg_C_diff),
# cor(micro_hydric$mass_g, micro_hydric$prop_above),
# body condition
cor(max_temps_hydric$SMI,
max_temps_hydric$chunk_avg_Tb_percentile_80_90),
cor(db_hydric$SMI, db_hydric$chunk_mean_Db_deg_C_diff),
cor(micro_hydric$SMI, micro_hydric$prop_above))) %>%
mutate(hydric_var = factor(hydric_var, levels = c("CEWL", "Plasma Osmolality",
#"Body Mass",
"Body Condition")))
corrs## hydric_var thermal_var R
## 1 CEWL Max Temp 0.03561470
## 2 CEWL Db 0.24284741
## 3 CEWL Micro Use -0.11107187
## 4 Plasma Osmolality Max Temp 0.16506018
## 5 Plasma Osmolality Db 0.03476896
## 6 Plasma Osmolality Micro Use -0.14174829
## 7 Body Condition Max Temp 0.26269876
## 8 Body Condition Db -0.21638149
## 9 Body Condition Micro Use 0.03188603Plot it:
removed body mass because body condition was more relevant
ggplot(corrs) +
aes(x = hydric_var,
y = thermal_var,
fill = R) +
geom_tile() +
scale_fill_distiller(palette = "RdYlBu", limits = c(-1, 1),
name = "Correlation\nCoefficient") +
savs_ggplot_theme2 +
scale_y_discrete(name = NULL, labels = c("Thermoregulatory\n Accuracy",
bquote("Maximum T"[b]),
"Microhabitat Use"
)) +
scale_x_discrete(name = NULL, labels = c("CEWL", "Plasma\nOsmolality",
#"Body\nMass",
"Body\nCondition")) +
# CEWL
annotate("text", x = 1, y = 1, label = "-0.03",
size = 2, family = "sans", color = "black") +
annotate("text", x = 1, y = 2, label = "-0.02",
size = 2, family = "sans", color = "black") +
annotate("text", x = 1, y = 3, label = "-0.01",
size = 2, family = "sans", color = "black") +
# plasma osmolality
annotate("text", x = 2, y = 1, label = "0.007*",
size = 3, family = "sans", color = "black") +
annotate("text", x = 2, y = 2, label = "0.01",
size = 2, family = "sans", color = "black") +
annotate("text", x = 2, y = 3, label = "-0.001",
size = 2, family = "sans", color = "black") +
# body mass
# annotate("text", x = 3, y = 1, label = "0.04",
# size = 2, family = "sans", color = "black") +
# annotate("text", x = 3, y = 2, label = "0.03",
# size = 2, family = "sans", color = "black") +
# annotate("text", x = 3, y = 3, label = "0.004",
# size = 2, family = "sans", color = "black") +
# body condition
annotate("text", x = 3, y = 1, label = "0.02",
size = 2, family = "sans", color = "black") +
annotate("text", x = 3, y = 2, label = "0.08*",
size = 3, family = "sans", color = "black") +
annotate("text", x = 3, y = 3, label = "0.005",
size = 2, family = "sans", color = "black") -> corr_plot
# arrange to get legend centered
ggarrange(corr_plot,
ncol = 1, nrow = 1,
common.legend = TRUE,
legend = "bottom"
) -> corr_plot_arrange
# export
ggsave(filename = "thermal_hydric_corr_plot.pdf",
plot = corr_plot_arrange,
path = "./results_figures",
device = "pdf",
dpi = 600,
units = "mm",
width = 100, height = 80)Hydric Conclusion
There are basically no detectable relationships between water balance (CEWL, osml, or mass) and thermoregulation or microhabitat use in BNLL.
Microhabitat Use ONLY
Prep Data
For logistic regression, I need to have the full, raw dataset (vs proportions by group for plots), so I will actually use the ‘microhab’ dataframe, created in the “Load Data” > “Microhabitat Use” section.
Basic Plots
ggplot(MH_use_across_indivs_sex) +
geom_bar(aes(x = date_chunks,
y = MH_use_proportion,
fill = microhabitat
),
position = "fill",
stat = "identity") +
theme_classic() +
facet_wrap(~sex_M_F)
ggplot(MH_use_across_indivs) +
geom_bar(aes(x = date_chunks,
y = MH_use_proportion,
fill = microhabitat
),
position = "fill",
stat = "identity") +
theme_classic()
Logistic Regression
Start with only ~ time period/chunk:
micro_glmm <- lme4::glmer(data = microhab,
above_below ~ date_chunks +
(1|individual_ID),
family="binomial")
summary(micro_glmm)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: above_below ~ date_chunks + (1 | individual_ID)
## Data: microhab
##
## AIC BIC logLik deviance df.resid
## 1688.7 1732.6 -836.4 1672.7 1767
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1617 -0.6023 -0.1103 0.4143 11.2158
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID (Intercept) 2.023 1.422
## Number of obs: 1775, groups: individual_ID, 39
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.709551 0.315963 -5.411 6.28e-08 ***
## date_chunksMay 10 - 20 0.419879 0.301382 1.393 0.164
## date_chunksMay 21 - 31 -0.004787 0.315245 -0.015 0.988
## date_chunksJun 1 - 11 0.429209 0.316201 1.357 0.175
## date_chunksJun 12 - 22 2.669814 0.307605 8.679 < 2e-16 ***
## date_chunksJun 23 - Jul 3 3.159663 0.291237 10.849 < 2e-16 ***
## date_chunksJul 4 - 13 4.177164 0.335053 12.467 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) d_M1-2 d_M2-3 d_J1-1 d_J1-2 d_2-J3
## dt_chM10-20 -0.484
## dt_chM21-31 -0.470 0.550
## dt_chnJ1-11 -0.477 0.537 0.560
## dt_chJ12-22 -0.516 0.550 0.565 0.642
## dt_chJ23-J3 -0.550 0.580 0.596 0.678 0.799
## dt_chnJ4-13 -0.490 0.505 0.519 0.592 0.715 0.784## contrast estimate SE df z.ratio p.value
## (Apr 26 - May 6) - (May 10 - 20) -0.41988 0.301 Inf -1.393 0.8058
## (Apr 26 - May 6) - (May 21 - 31) 0.00479 0.315 Inf 0.015 1.0000
## (Apr 26 - May 6) - (Jun 1 - 11) -0.42921 0.316 Inf -1.357 0.8244
## (Apr 26 - May 6) - (Jun 12 - 22) -2.66981 0.308 Inf -8.679 <.0001
## (Apr 26 - May 6) - (Jun 23 - Jul 3) -3.15966 0.291 Inf -10.849 <.0001
## (Apr 26 - May 6) - (Jul 4 - 13) -4.17716 0.335 Inf -12.467 <.0001
## (May 10 - 20) - (May 21 - 31) 0.42467 0.293 Inf 1.450 0.7743
## (May 10 - 20) - (Jun 1 - 11) -0.00933 0.298 Inf -0.031 1.0000
## (May 10 - 20) - (Jun 12 - 22) -2.24993 0.289 Inf -7.784 <.0001
## (May 10 - 20) - (Jun 23 - Jul 3) -2.73978 0.272 Inf -10.081 <.0001
## (May 10 - 20) - (Jul 4 - 13) -3.75728 0.318 Inf -11.813 <.0001
## (May 21 - 31) - (Jun 1 - 11) -0.43400 0.296 Inf -1.465 0.7658
## (May 21 - 31) - (Jun 12 - 22) -2.67460 0.291 Inf -9.206 <.0001
## (May 21 - 31) - (Jun 23 - Jul 3) -3.16445 0.273 Inf -11.577 <.0001
## (May 21 - 31) - (Jul 4 - 13) -4.18195 0.319 Inf -13.090 <.0001
## (Jun 1 - 11) - (Jun 12 - 22) -2.24060 0.264 Inf -8.488 <.0001
## (Jun 1 - 11) - (Jun 23 - Jul 3) -2.73045 0.245 Inf -11.156 <.0001
## (Jun 1 - 11) - (Jul 4 - 13) -3.74796 0.295 Inf -12.723 <.0001
## (Jun 12 - 22) - (Jun 23 - Jul 3) -0.48985 0.191 Inf -2.569 0.1356
## (Jun 12 - 22) - (Jul 4 - 13) -1.50735 0.244 Inf -6.184 <.0001
## (Jun 23 - Jul 3) - (Jul 4 - 13) -1.01750 0.210 Inf -4.848 <.0001
##
## Results are given on the log odds ratio (not the response) scale.
## P value adjustment: tukey method for comparing a family of 7 estimatespairwise comparison groups for when I make figures:
- Apr 26 - May 6 = A
- May 10 - 20 = A
- May 21 - 31 = A
- Jun 1 - 11 = A
- Jun 12 - 22 = B
- Jun 23 - Jul 3 = B
- Jul 4 - 13 = C (wow!!)
Compare a model with only ~ chunks vs chunks*sex:
micro_glmm_sex <- lme4::glmer(data = microhab,
above_below ~ date_chunks*sex_M_F +
(1|individual_ID),
family="binomial")
summary(micro_glmm_sex)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: above_below ~ date_chunks * sex_M_F + (1 | individual_ID)
## Data: microhab
##
## AIC BIC logLik deviance df.resid
## 1677.1 1759.4 -823.6 1647.1 1760
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0934 -0.5633 -0.0851 0.4084 14.0614
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID (Intercept) 1.828 1.352
## Number of obs: 1775, groups: individual_ID, 39
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.4451 0.4263 -3.390 0.000699 ***
## date_chunksMay 10 - 20 1.1237 0.4159 2.702 0.006899 **
## date_chunksMay 21 - 31 0.6417 0.4552 1.410 0.158613
## date_chunksJun 1 - 11 1.3148 0.4952 2.655 0.007936 **
## date_chunksJun 12 - 22 2.8696 0.4985 5.756 8.61e-09 ***
## date_chunksJun 23 - Jul 3 3.4974 0.4714 7.420 1.17e-13 ***
## date_chunksJul 4 - 13 3.4882 0.5530 6.307 2.84e-10 ***
## sex_M_FMale -0.4569 0.6125 -0.746 0.455696
## date_chunksMay 10 - 20:sex_M_FMale -1.5098 0.6171 -2.447 0.014416 *
## date_chunksMay 21 - 31:sex_M_FMale -1.2489 0.6350 -1.967 0.049199 *
## date_chunksJun 1 - 11:sex_M_FMale -1.5622 0.6467 -2.416 0.015706 *
## date_chunksJun 12 - 22:sex_M_FMale -0.4402 0.6282 -0.701 0.483534
## date_chunksJun 23 - Jul 3:sex_M_FMale -0.6205 0.5934 -1.046 0.295746
## date_chunksJul 4 - 13:sex_M_FMale 0.6935 0.6864 1.010 0.312392
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.0650045 (tol = 0.002, component 1)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: above_below
## Chisq Df Pr(>Chisq)
## date_chunks 283.6728 6 < 2.2e-16 ***
## sex_M_F 5.3173 1 0.021115 *
## date_chunks:sex_M_F 20.3379 6 0.002411 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: microhab
## Models:
## micro_glmm: above_below ~ date_chunks + (1 | individual_ID)
## micro_glmm_sex: above_below ~ date_chunks * sex_M_F + (1 | individual_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## micro_glmm 8 1688.7 1732.6 -836.36 1672.7
## micro_glmm_sex 15 1677.2 1759.4 -823.57 1647.2 25.581 7 0.0005982 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Okay, adding sex was technically a better model based on some metrics, but how many time periods was microhabitat use actually different between the sexes?
emmeans_micro_sex_glmm_sex_diffs <- data.frame(emmeans(micro_glmm_sex,
pairwise ~ sex_M_F | date_chunks)$contrasts) %>%
dplyr::select(date_chunks, contrast, p.value)
emmeans_micro_sex_glmm_sex_diffs## date_chunks contrast p.value
## 1 Apr 26 - May 6 Female - Male 0.455695676
## 2 May 10 - 20 Female - Male 0.001642991
## 3 May 21 - 31 Female - Male 0.007776951
## 4 Jun 1 - 11 Female - Male 0.001975646
## 5 Jun 12 - 22 Female - Male 0.150030024
## 6 Jun 23 - Jul 3 Female - Male 0.066251295
## 7 Jul 4 - 13 Female - Male 0.726818004May - June, microhabitat use was significantly different based on sex. Females were spending more time in burrows than males (based on the plots above).
To keep it simple, let’s keep the main figure to showing both males and females aggregated over time, but I can also make a secondary figure showing their differences in microhabitat use specifically in May-early June.
Probabilities from GLMM
Now, I want to use the model estimates/coefficients to calculate the predicted probability of lizards being belowground. First, save the estimates and standard errors. I can do this by getting the model-estimated means. Then, I can use this equation to calculate the model-predicted probability of being belowground:
p = exp(emmean) / (1 + exp(emmean))
# one set for M-F differences
emmeans_probabilities_MF <- data.frame(emmeans(micro_glmm_sex,
pairwise ~ sex_M_F | date_chunks)$emmeans) %>%
# calculate prob
mutate(exp_emmean = exp(emmean),
prob = exp_emmean / (1 + exp_emmean),
# and confidence intervals, in case helpful
exp_LCL = exp(asymp.LCL),
prob_LCL = exp_LCL / (1 + exp_LCL),
exp_UCL = exp(asymp.UCL),
prob_UCL = exp_UCL / (1 + exp_UCL)
) %>%
# only save the probabilities for when M & F were sig diff
dplyr::filter(date_chunks %in% c("May 10 - 20", "May 21 - 31", "Jun 1 - 11"))
emmeans_probabilities_MF## sex_M_F date_chunks emmean SE df asymp.LCL asymp.UCL exp_emmean
## 1 Female May 10 - 20 -0.3214039 0.4264786 Inf -1.157286 0.51447873 0.72513033
## 2 Male May 10 - 20 -2.2880805 0.4562992 Inf -3.182410 -1.39375044 0.10146103
## 3 Female May 21 - 31 -0.8034043 0.4608742 Inf -1.706701 0.09989252 0.44780191
## 4 Male May 21 - 31 -2.5091531 0.4450839 Inf -3.381502 -1.63680471 0.08133709
## 5 Female Jun 1 - 11 -0.1303629 0.4964374 Inf -1.103362 0.84263653 0.87777684
## 6 Male Jun 1 - 11 -2.1494058 0.4224566 Inf -2.977406 -1.32140601 0.11655339
## prob exp_LCL prob_LCL exp_UCL prob_UCL
## 1 0.42033365 0.31433798 0.23916069 1.6727663 0.6258558
## 2 0.09211495 0.04148553 0.03983304 0.2481429 0.1988097
## 3 0.30929778 0.18146343 0.15359208 1.1050521 0.5249524
## 4 0.07521900 0.03399637 0.03287862 0.1946009 0.1629003
## 5 0.46745536 0.33175375 0.24911043 2.3224822 0.6990202
## 6 0.10438676 0.05092478 0.04845712 0.2667600 0.2105845# one set for overall probabilities
emmeans_probabilities <- data.frame(emmeans(micro_glmm, "date_chunks")) %>%
# calculate prob
mutate(exp_emmean = exp(emmean),
prob = exp_emmean / (1 + exp_emmean),
# and confidence intervals, in case helpful
exp_LCL = exp(asymp.LCL),
prob_LCL = exp_LCL / (1 + exp_LCL),
exp_UCL = exp(asymp.UCL),
prob_UCL = exp_UCL / (1 + exp_UCL)
) %>%
# only save the probabilities for when M & F were the same
dplyr::filter(date_chunks %nin% c("May 10 - 20", "May 21 - 31", "Jun 1 - 11")
) %>%
# add "sex" col
mutate(sex_M_F = "Both")
emmeans_probabilities## date_chunks emmean SE df asymp.LCL asymp.UCL exp_emmean
## 1 Apr 26 - May 6 -1.7095505 0.3159632 Inf -2.3288271 -1.090274 0.1809471
## 2 Jun 12 - 22 0.9602636 0.3067858 Inf 0.3589745 1.561553 2.6123850
## 3 Jun 23 - Jul 3 1.4501123 0.2887028 Inf 0.8842652 2.015959 4.2635933
## 4 Jul 4 - 13 2.4676138 0.3292154 Inf 1.8223634 3.112864 11.7942692
## prob exp_LCL prob_LCL exp_UCL prob_UCL sex_M_F
## 1 0.1532220 0.09740994 0.08876349 0.3361244 0.2515667 Both
## 2 0.7231746 1.43186024 0.58879216 4.7662161 0.8265760 Both
## 3 0.8100157 2.42120453 0.70770529 7.5079274 0.8824626 Both
## 4 0.9218400 6.18646237 0.86084948 22.4853524 0.9574203 Both# put them together
barplot_probs <- emmeans_probabilities %>%
rbind(emmeans_probabilities_MF) %>%
dplyr::select(date_chunks, sex_M_F,
prob, prob_LCL, prob_UCL) %>%
mutate(sex_M_F = factor(sex_M_F))
summary(barplot_probs)## date_chunks sex_M_F prob prob_LCL
## Apr 26 - May 6:1 Both :4 Min. :0.07522 Min. :0.03288
## May 10 - 20 :2 Female:3 1st Qu.:0.11660 1st Qu.:0.05853
## May 21 - 31 :2 Male :3 Median :0.36482 Median :0.19638
## Jun 1 - 11 :2 Mean :0.40771 Mean :0.30091
## Jun 12 - 22 :1 3rd Qu.:0.65924 3rd Qu.:0.50387
## Jun 23 - Jul 3:1 Max. :0.92184 Max. :0.86085
## Jul 4 - 13 :1
## prob_UCL
## Min. :0.1629
## 1st Qu.:0.2208
## Median :0.5754
## Mean :0.5340
## 3rd Qu.:0.7947
## Max. :0.9574
## Pretty Barplot
ggplot() +
# raw observed proportions
geom_bar(data = MH_use_across_indivs,
aes(x = date_chunks,
y = MH_use_proportion,
fill = microhabitat
),
position = "fill",
stat = "identity"
) +
# model-predicted probabilities
geom_errorbar(data = barplot_probs,
aes(
x = date_chunks,
ymin = prob_LCL,
ymax = prob_UCL,
),
color = "black",
width = 0.1
) +
geom_point(data = barplot_probs,
aes(
x = date_chunks,
y = prob,
shape = sex_M_F,
),
color = "black",
size = 2
) +
# labels
xlab("Time Period") +
ylab("Microhabitat Use (proportion of time)") +
# hydric phys msmsts
geom_vline(xintercept = 0.5, # first weekend of msmts
size = 0.5, # default size is 0.5
linetype = "dashed") +
geom_vline(xintercept = 1.5, # second weekend of msmts
size = 0.5,
linetype = "dashed") +
geom_vline(xintercept = 7.5, # last few msmts
size = 0.5,
linetype = "dashed") +
# pretties
scale_shape_manual(values = c(17, 16, 15), name = "Lizard Sex") +
scale_fill_manual(values = my_brew_cols_4, name = "Microhabitat",
labels = c("Open", "Partial Shade", "Full Shade", "Burrow")) +
scale_x_discrete(labels = c("Apr 26\n-May 6",
"May\n10-20",
"May\n21-31",
"Jun\n1-11",
"Jun\n12-22",
"Jun 23\n-Jul 3",
"Jul\n4-13"
),
name = NULL) +
scale_y_continuous(limits = c(0, 1.1), expand = c(0,0),
breaks = c(seq(0, 1, by = 0.2)),
labels = c(seq(0, 1, by = 0.2))) +
# significance annotation for time
annotate(geom = "text", x = 1, y = 1.08, label = "a", size = 3) +
annotate(geom = "text", x = 2, y = 1.08, label = "a", size = 3) +
annotate(geom = "text", x = 3, y = 1.08, label = "a", size = 3) +
annotate(geom = "text", x = 4, y = 1.08, label = "a", size = 3) +
annotate(geom = "text", x = 5, y = 1.08, label = "b", size = 3) +
annotate(geom = "text", x = 6, y = 1.08, label = "b", size = 3) +
annotate(geom = "text", x = 7, y = 1.08, label = "c", size = 3) +
# significance annotation for sex
annotate(geom = "text", x = 2, y = 1.015, label = "*", size = 5) +
annotate(geom = "text", x = 3, y = 1.015, label = "*", size = 5) +
annotate(geom = "text", x = 4, y = 1.015, label = "*", size = 5) +
savs_ggplot_theme -> fancy_micro_plot
fancy_micro_plot
# arrange to get legend centered
ggarrange(fancy_micro_plot,
ncol = 1, nrow = 1,
common.legend = TRUE,
legend = "right"
) -> fancy_micro_plot_arrange
ggsave(filename = "microhabitat_use.pdf",
plot = fancy_micro_plot_arrange,
path = "./results_figures",
device = "pdf",
dpi = 600,
units = "mm",
width = 160, height = 80)